Vortex proliferation in the Berezinskii-Kosterlitz-Thouless regime on a 
two-dimensional lattice of Bose-Einstein condensates 
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We observe the proliferation of vortices in the Berezinskii-Kosterlitz-Thouless regime on a two- 
dimensional array of Josephson-coupled Bose-Einstein condensates. As long as the Josephson (tun- 
neling) energy J exceeds the thermal energy T, the array is vortex- free. With decreasing J/T, 
vortices appear in the system in ever greater numbers. We confirm thermal activation as the vortex 
formation mechanism and obtain information on the size of bound vortex pairs as J/T is varied. 
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One of the defining characteristics of superfluids is 
long-range phase coherence [1], which may be destroyed 
by quantum fluctuations, as in the Mott-insulator tran- 
sition @, 0], or thermal fluctuations, e.g. in one- 
dimensional Bose gases [4, 5] and in a double- well system 
[f|. In two dimensions (2D), Berezinskii Kosterlitz 
and Thouless [8] (BKT) developed an elegant description 
of thermal phase fluctuations based on the unbinding of 
vortex-antivortex pairs, i.e. pairs of vortices of opposite 
circulation. The BKT picture applies to a wide variety 
of 2D systems, among them Josephson junction arrays 
(JJA), i.e. arrays of superfluids in which phase coherence 
is mediated via a tunnel coupling J between adjacent 
sites. Placing an isolated (free) vortex into a JJA is ther- 
modynamically favored if its free energy F = E—TS < 0. 
In an array of period d the vortex energy diverges with ar- 
ray size R as E w J\og(R/d) 0], but may be offset by an 
entropy gain S w \og(R/d) due to the available w R 2 /d 2 
sites. This leads to a critical condition (J/T) cr n ~ 1 in- 
dependent of system size, below which free vortices will 
proliferate. In contrast, tightly bound vortex-antivortex 
pairs are less energetically costly and show up even above 
(«7/T) cr j{. The overall vortex density is thus expected to 
grow smoothly with decreasing J/T in the BKT crossover 
regime. 

Transport measurements, both in continuous superflu- 
ids [13, HH and superconducting JJA [l2j have confirmed 
the predictions of BKT, without however directly observ- 
ing its microscopic mechanism, vortex-antivortex unbind- 
ing. A recent beautiful experiment [131 ] in a continuous 
2D Bose gas measured the phase-phase decay function 
through the BKT cross-over, and saw evidence for ther- 
mal vortex formation. For related theoretical studies see 
e.g. jlij . In this work we present more detailed vortex- 
formation data, collected in a 2D array of BECs with ex- 
perimentally controllable Joephson couplings. The sys- 
tem was studied theoretically in [l5l ]. 

Our experiment starts with production of a partially 
Bose-condensed sample of 87 Rb atoms in a harmonic, ax- 
ially symmetric magnetic trap with oscillation frequen- 
cies {w Xj y,u) z } = 27r{6.95, i5.0}Hz. The number of con- 
densed atoms is kept fixed around 6 x I0 5 as the tem- 
perature is varied. We then transform this system into 
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FIG. 1: Experimental system, (a) 2D optical lattice intensity 
profile. A lattice of Josephson-coupled BECs is created in the 
white-shaded area. The central box marks the basic building 
block of our system, the double- well potential shown in (b). 
The barrier height Vol and the number of condensed atoms 
per well, N we u , control the Josephson coupling J, which acts 
to lock the relative phase A(j>. A cloud of uncondensed atoms 
at temperature T induces thermal fluctuations and phase de- 
fects in the array when J < T. (c) Experimental sequence: A 
BEC (i) is loaded into the optical lattice over 10 s, suppress- 
ing J to values around T. We allow 2 s for thermalization. To 
probe the system, we ramp off the lattice on a faster timescale 
t r [23| and take images of the recombined condensate. When 
J is reduced below T (ii)-(vii), vortices (dark spots) appear 
as remnants of the thermal fluctuations in the array. 



a Josephson junction array, as illustrated in Fig. [T] In 
a 10 s linear ramp, we raise the intensity of a 2D hexag- 
onal optical lattice [16[ of period d — 4.7/im in the x-y 
plane. The resulting potential barriers of height Vol be- 
tween adjacent sites [FigQJb)] rise above the condensate's 
chemical potential around Vol ~ 250 — 300 Hz, splitting- 
it into an array of condensates which now communicate 



only through tunneling. This procedure is adiabatic even 
with respect to the longest-wavelength phonon modes of 
the array [H,[l^] over the full range of Vol in our exper- 
iments. Each of the sw 190 occupied sites (15 sites across 
the BEC diameter 2 x Rtf ~ 68 fim [20() now contains a 
macroscopic BEC, with N we ii w 7000 condensed atoms in 
each of the central wells at a temperature T that can be 
adjusted between 30— 70 nK. By varying Vol in a range 
between 500 Hz and 2 kHz we tune J between 1.5 [iK 
and 5nK, whereas the "charging" energy E c , defined in 
[l|, is on the order of a few pK, much smaller than both 
J and T. In this regime, thermal fluctuations of the rela- 
tive phases A<j)Th ~ \fT/ J are expected, while quantum 
fluctuations A</>q w (Ec/AJ) 1 ' 4 are neglig ible Q. 

The suppression of the Josephson coupling greatly sup- 
presses the energy cost of phase fluctuations in the x-y 
plane, between condensates, J[l — cos(A<^)], compared to 
the cost of axial (z) phase fluctuations inside the con- 
densates [2l|. As a result, axial phase fluctuations re- 
main relatively small, and each condensate can be ap- 
proximated as a single-phase object [22| . 

After allowing 2 s for thermalization, we initiate our 
probe sequence. We first take a nondestructive thermom- 
etry image in the x-z plane, from which the temperature 
T and, from the axial condensate size R z , the number 
of condensed particles per well, N we u, is obtained (see 
below). To observe the phase fluctuations we then turn 
down the optical lattice on a time-scale t r [23|, which 
is fast enough to trap phase winding defects, but slow 
enough to allow neighboring condensates to merge, pro- 
vided their phase difference is small. Phase fluctuations 
are thus converted to vortices in the reconnected conden- 
sate, as has been observed in the experiments of Scherer 
et al. (24|. We then expand the condensate by a factor 
of 6 and take a destructive image in the x-y plane. 

Figure [TJc) illustrates our observations: (ii)-(vii) is a 
sequence of images at successively smaller J/T (measured 
in the center of the array (25[). Vortices, with their cores 
visible as dark "spots" in (iii)-(vii), occur in the BEC 
center around J/T = 1. Vortices at the BEC edge ap- 
pear earlier, as here the magnetic trap potential adds 
to the tunnel barrier, suppressing the local J/T below 
the quoted value. That the observed "spots" are indeed 
circulation-carrying vortices and antivortices is inferred 
from their slow w 100 ms decay after the optical lattice 
ramp-down, presumably dominated by vortex-antivortex 
annihilation. From extensive experiments on vortices in 
our system we know that circulation-free "holes" fill so 
quickly due to positive mean field pressure, that they 
do not survive the pre-imaging expansion. Vortices with 
identical circulation would decay by dissipative motion 
to the BEC edge, in our trap over > 10 s. 

To investigate the thermal nature of phase fluctua- 
tions, we study vortex activation while varying J at 
different temperatures. For a quantitative study, accu- 
rate parameter estimates are required. The Josephson- 
coupling energy J is obtained from 3D numerical sim- 
ulations of the Gross-Pitaevskii equation (CPE) for 
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FIG. 2: Quantitative study of vortices. The areal density of 
vortices is quantified by the plotted T> defined in the text. T> is 
extracted only from the central 11% of the condensate region 
[circle in inset (a)] to minimize effects of spatial inhomogene- 
ity. (b) T> vs J for two datasets with distinct "cold" and "hot" 
temperatures. Each point represents one experimental cycle. 
The increase in T> with decreasing J < 100 nK signals the 
spontaneous appearance of vortices, while the "background" 
T> < 0.01 for J > 200 nK is not associated with vortices. 
Vortices clearly proliferate at larger J for the "hot" data, 
indicating thermal activation as the underlying mechanism. 
The large scatter in D at low J is due to shot noise on the 
small average number of vortices in the central condensate 
region, (c) same data as in (b), but averaged within bins of 
size A[log(J)] = 0.15. Error bars of T> are standard errors, 
(d) same data as (b), but plotted vs J/T. "Cold" and "hot" 
datasets almost overlap on what appears to be a universal 
vortex activation curve, as confirmed by averaging [inset (e)], 
clearly revealing the underlying competition of J and T. 



the central double- well system [y, [2(| [FigQJb)], se lf- 
consistently including mean-field interactions of both 
condensed and uncondensed atoms (27j. A use- 
ful approximation for J in our experiments is [25| : 
J(V L,N weU ,T) » N weU x 0.315n#exp[A„ e „/3950 - 
V O L/2UHz](l + 0.59T/100nK). The finite-T correction 
to J arises from both the lifting-up of the BEC's chem- 
ical potential and the axial compression by the thermal 
cloud's repulsive mean field, but does not take into ac- 
count the effects of phase fluctuations on J (in condensed- 
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matter language, we calculate the bare J). N we u is de- 
termined by comparison of the experimentally measured 
R z , to R z (VoL,N we ii,T) obtained from GPE simula- 
tions. Both experimental and simulated R z are obtained 
from a fit to the distribution of condensed and uncon- 
densed atoms, to a Thomas-Fermi profile plus mean-field- 
modified Bose function 27]. In determination of all J 
values, there is an overall systematic multiplicative un- 
certainty A J I J = *1.6, dominated by uncertainties in 
the optical lattice modulation contrast, the absolute in- 
tensity calibration, and magnification in the image used 
to determine N we u. In comparing J for "hot" and "cold" 
clouds (see Fig. [2]) there is a relative systematic error of 
15% associated with image fitting and theory uncertain- 
ties in the thermal-cloud mean-field correction to J. 

The qualitative results of our work are consistent 
whether we use an automated vortex-counting routine 
or count vortices by hand, but the former shows signs of 
saturation error at high vortex density, and the latter is 
vulnerable to subjective bias. As a robust vortex-density 
surrogate we therefore use the "roughness" T> of the con- 
densate image caused by the vortex cores. Precisely, we 
define V as the normalized variance of the measured col- 
umn density profile from a fit to a smooth finite-T Bose 
profile [273, with a small constant offset subtracted to 
account e.g. for imaging noise. To limit spatial inhomo- 
geneity in J, caused by spatially varying condensate den- 
sity and optical lattice intensity, to < 10%, V is extracted 
only from the central 11% of the condensate area which 
contains 20 lattice sites [Fig. Ufa)]. Comparison to au- 
tomated vortex-counts shows that T> is roughly linear in 
the observed number of vortices, irrespective of the sign 
of their circulation, with a sensitivity of ss 0.01/vortex. 

Figure [5] shows results of our quantitative study. In 
Fig. [21b), we plot T> vs J for two datasets with dis- 
tinct temperatures. At large J > 200 nK a background 
T> < 0.01 is observed, that is not associated with vortices, 
but due to residual density ripples remaining after the 
optical lattice ramp-down. Vortex proliferation, signaled 
by a rise of T> above 0.01, occurs around J s» 100 nK 
for "hot" BECs and at a distinctly lower J ss 50 nK for 
"cold" BECs [confirmed by the averaged data shown in 
Fig. [UJc)], indicating thermal activation as the vortex 
formation mechanism. Plotting the same data vs J/T 
in Fig. HJd) shows collapse onto a universal vortex acti- 
vation curve, providing strong evidence for thermal ac- 
tivation. A slight residual difference becomes visible in 
the averaged "cold" vs "hot" data [FigJSJe)] , perhaps be- 
cause of systematic differences in our determination of J 
at different temperatures. 

The vortex density T> by itself provides no distinction 
between bound vortex-antivortex pairs and free vortices. 
In the following we exploit the flexibility of optical po- 
tentials to distinguish free or loosely bound vortices from 
tightly bound vortex-antivortex pairs. A "slow" optical 
lattice ramp-down allows time for tightly bound pairs 
to annihilate before they can be imaged. By slowing 
down the ramp-down duration r [inset of Fig. [3] (a)], 
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FIG. 3: (a) Vortex density T> probed at different optical lat- 
tice ramp-down timescales r. A slow ramp provides time 
for tightly bound vortex-antivortex pairs to annihilate, allow- 
ing selective counting of loosely bound or free vortices only, 
whereas a fast ramp probes both free and tightly bound vor- 
tices. A fit to the vortex activation curve determines its mid- 
point (J/T) 50% , its 27% - 73% width A(J/T) 27 -73, and the 
limiting values T> < (£>>) well below (above) (J/T) 50 %. (b) 
A downshift in (J/T) 50% is seen for slow ramps, consistent 
with the occurrence of loosely bound or free vortices at lower 
J/T only, (c) Mapping between ramp-down timescale r and 
estimated size of the smallest pairs surviving the ramp (up- 
per axis). The difference T> < — T) > measures the number of 
vortices surviving the ramp (right axis). Comparison to sim- 
ulated vortex distributions yields a size estimate of the small- 
est surviving pairs (upper axis). Inset: smallest possible pair 
sizes in a hexagonal array, I: d/\/3, II: d, III: Id/yfZ. 



we therefore selectively probe vortices at increasing spa- 
tial scales. This represents an attempt to approach the 
"true" BKT vortex unbinding crossover that is comple- 
mentary to transport measurements employed so success- 
fully in superconductive and liquid Helium systems. 

Figure [3](a) shows vortex activation curves, probed 
with two different ramp-down times. Two points are 
worth noticing: First, a slow ramp compared to a fast 
one shows a reduction of the vortex density T> < in arrays 
with fully randomized phases at low J/T. The difference 
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directly shows the fraction of tightly bound pairs that an- 
nihilate on the long ramp. Second, a slower ramp shows 
vortex activation at lower (J/T) 50 %, confirming that free 
or very loosely bound vortices occur only at higher T 
(lower J). Specifically, the data clearly show a range 
around J/T « 1.4 where only tightly bound pairs exist. 
Figure [3fb) quantitatively shows the shift of (J/T) 50% 
from 1.4 to 1.0 with slower ramp time. We can make 
a crude mapping of the experimental ramp-down time- 
scale to theoretically more accessible vortex-antivortex 
pair sizes as follows: In Fig. [3jc), we see the decrease of 
the saturated (low- J/T) vortex density T> < with increas- 
ing ramp timescale r. The right axis shows the inferred 
number of vortices that survived the ramp. We compare 
this number of surviving vortices to simulations [28| of 
a 20-site hexagonal array with random phases. In these 
simulations we find, on average, a total of 10 vortices, 
6 of which occur in nearest-neighbor vortex-antivortex 
pairs [configuration I in Fig. EJc)], 1.7 (0.4) occur in 
configuration II (III) respectively, and 1.9 occur in larger 
pairs or as free vortices. Experimentally 11 vortices are 
observed for the fastest ramps, in good agreement with 
the expected total number of vortices. For just some- 
what slower ramps of t w 5 ms, only 3 vortices survive, 
consistent with only vortices in configuration II & III or 
larger remaining (indicated in Fig. 02 top axis) [29j |. For 
t > 30 ms ramps less than 2 vortices remain, according 
to our simulations spaced by more than 2<i/v / 3- Thus we 
infer that ramps of r ~ 30 ms or longer allow time for 
bound pairs of spacing < 2c?/ v3 to decay before we ob- 



serve them. The downward shift of (J/T) 50% in Fig. \3j[b) 
thus tells us that loosely bound pairs of size larger than 
2d/V3, or indeed free vortices, do not appear in quan- 
tity until J/T < 1.0, whereas more tightly bound vortex 
pairs appear in large number already for J/T < 1.4. 

A further interesting observation concerns the width 
of the vortex activation curve. The relative width, de- 
termined from fits to data such as the ones shown in 
Fig. Ha), is A(J/T) 27 _ 73 /(J/T) 50% » 0.3, independent 
of ramp-down duration. This width is neither as broad 
as in a double- well system @, |3(J, where the coherence 
factor rises over a range A( J/T) 2 7_73/( J/T) 50 % « 1.4, 
nor as broad as expected from our simulations [28| of 
an array of uncoupled phase s, ea ch fluctuating inde- 
pendently with A(f>RMS = \/T I J, for which we find 
A(J/T) 2 7-73/{J/T) 5Q% w 0.85. Presumably collective 
effects in the highly multiply connected lattice narrow 
the curve. On the other hand, the width is 3 times larger 
than the limit due to spatial inhomogeneity in J, suggest- 
ing contributions to the width due to finite-size effects 
or perhaps revealing the intrinsically smooth behavior of 
vortex activation in the BKT regime. 

In conclusion, we have probed vortex proliferation in 
the BKT regime on a 2D lattice of Josephson-coupled 
BECs. Allowing variable time for vortex-antivortex pair 
annihilation before probing the system provides a time- 
to-length mapping, which reveals information on the size 
of pairs with varying J/T. We acknowledge illuminating 
conversations with Leo Radzihovsky and Victor Gurarie. 
This work was funded by NSF and NIST. 



[*] Quantum Physics Division, National Institute of Stan- 
dards and Technology. 

[1] A. Leggett, Rev. Mod. Phys. 73, 307 (2001). 

[2] M. Greiner et al, Nature 415, 39 (2002). 

[3] D. Jaksch et al, Phys. Rev. Lett. 81, 003108 (1998). 

[4] S. Richard et al, Phys. Rev. Lett. 91, 010405 (2003); 

[5] D. Hellweg et al, Phys. Rev. Lett. 91, 010406 (2003). 

[6] R. Gati et al, Phys. Rev. Lett. 96, 130404 (2006). 

[7] V. Berezinskii, Sov. Phys.-JETP 32, 493 (1971); 34, 610 
(1972). 

[8] J. Kosterlitz, D. Thouless, J. Phys. C 6, 1181 (1973). 
[9] M. Tinkham, Introduction to Superconductivity, 
McGraw-Hill, Inc., New York (1996). 
[10] G. Agnolet et al, Phys. Rev. B 39, 8934 (1989). 
[11] A. Safonov et al, Phys. Rev. Lett. 81, 4545 (1998) 
[12] D. Resnick et al, Phys. Rev. Lett. 47, 1542 (1981). 
[13] Z. Hadzibabic et al, Nature 441, 1118 (2006). 
[14] T. Simula and P. Blakie, Phys. Rev. Lett. 96, 020404 
(2006). 

[15] A. Trombettoni et al, New J. Phys. 7, 57 (2005). 

[16] Three circularly polarized laser beams (A = 810. lnm) in- 
tersect in a tripodlike configuration, with 6 — 6.6° angles 
to the z-axis. Calculation of the optical dipole potential 
[d| includes counterrotating terms and interaction with 
both the Dl and D2 lines, as well as the 'fictitious mag- 
netic field' due to the circular polarization. The tilted 
bias field of the TOP trap makes V « 0.5[I3]- 



[17] R. Grimm et al, Adv. At. Mol. Opt. Phys. 42, 95 (2000). 

[18] J. Javanainen, Phys. Rev. A 60, 4902 (1999). 

[19] K. Burnett et al, J. Phys. B 35, 1671 (2002). 

[20] To avoid radial flows during Vol ramp-up, Rtf is kept 
constant by balancing the lattice-enhanced mean field 
pressure with radial confinement due to the optical lattice 
envelope, by the choice of a 67pm 1/e 2 intensity waist. 
Axial (z) confinement is due to the magnetic trap alone. 

[21] D. Petrov et al, Phys. Rev. Lett. 87, 050404 (2001). 

[22] In the axial condensate region between z — 
— Rz/3, +R z /3, where according to our 3D GPE 
simulations 85% of the tunnel current is localized 
and hence the relative phase is measured, axial phase 
fluctuations [21[ vary between « 600 mrad ("cold" data 
in Fig. and 800 mrad ("hot") in the regime J/T « 1. 

[23] Within a dataset, the ramp-down rate is kept fixed, t r = 
t x Vol/1-3 kH z, t = 18 ms if not otherwise indicated. 

[24] D. Scherer et al, Phys. Rev. Lett. 98, 110402 (2007). 

[25] J is averaged over junctions within the central 11% of 
the array area, from which all quantitative experimental 
results are extracted. 

[26] D. Ananikian and T. Bergeman, Phys. Rev. A 73, 013604 
(2006). 

[27] M. Naraschewski and D. Stamper-Kurn, Phys. Rev. A 

58, 2423 (1998). 
[28] Following [2^], in simulations we count a vortex if all 

three phase differences in an elemental triangle of junc- 



tions are £ (0, 71"), or if all are 6 (— 7r,0). core in the bulk condensate (after Vol rampdown). 

The very short annihilation time of configuration I pairs [30] L. Pitaevskii and S. Stringari, Phys. Rev. Lett. 

[ r < 5ms in Fig.[3jc) ] is not unexpected: their spacing 180402 (2001). 
d/y/3 = 2.8 fxm is comparable to the diameter of a vortex 



